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Abstract. We analytically solve two problems that may be useful in the context of the recent observation of 
matter wave bright solitons in a one-dimensional attractive atomic Bose gas. The first problem is strictly 
beyond mean field: From the Bethe ansatz solution we extract the internal correlation function of the 
particle positions in the quantum soliton, that is for a fixed center of mass position. The second problem 
is solved in the limit of a large number of particles, where the mean field theory is asymptotically correct: 
It deals with the number of excitations created by the opening of the trap, starting from a pure soliton in 
a weakly curved harmonic potential. 

PACS. 03.75.Lm Tunneling, Josephson effect, Bose-Einstein condensates in periodic potentials, solitons, 
vortices, and topological excitations - 03.75.Hh Static properties of condensates; thermodynamical, statis- 
tical, and structural properties - 03.75.Kk Dynamic properties of condensates; collective and hydrodynamic 
excitations, superfluid flow 



Experiments with cold atoms have now acquired a high 
degree of control of the key parameters of the system. Us- 
ing transverse confinement of the atoms by non-dissipative 
optical potentials it is possible to freeze the atomic mo- 
tion along one or several directions, realizing in this way 
quantum gases with reduced dimensionality pQ. Further- 
more, thanks to Feshbach resonances driven by a magnetic 
field, one can adjust almost at will the interaction strength 
between the atoms [2] . The combination of these two ex- 
perimental tools has recently allowed the observation of 
bright solitons in a one-dimensional (ID) Bose gas, either 
a single soliton [3] or a train of solitons [J] . 

This leads to a renewed interest [51El[7j in the ID non- 
relativistic Bose gas on the free line with zero range at- 
tractive interactions, that is with a binary interaction po- 
tential modeled by a Dirac delta with a negative coupling 
constant g, V{x\ — X2) = gS(xi — X2), a model that is an 
acceptable approximation of reality under conditions de- 
fined in [5] . On a theoretical point of view, it is well known 
that eigenstates and eigenenergies of the corresponding TV- 
body Hamiltonian may be obtained by the Bethe ansatz: 
Historically the focus was put mainly on the study of the 
ground state wavefunction <f)(xi, . . . ,xn) [WlO]. which is 
a collective bound state of the N particles, the so-called 
iV-particle quantum soliton, with a delocalized center of 
mass of vanishing momentum. The extension of the Bethe 
ansatz to excited states is however possible, and one finds 
that the generic excited state corresponds to a set of quan- 
tum solitons with arbitrary atom numbers and different 
momenta per particle [TT] . A key property to keep in mind 
is the full separability of the center of mass variables and 



the internal variables of the gas (e.g. the relative coordi- 
nates of the particles) , which holds since the gas is on the 
free line with open boundary conditions. 

Despite the knowledge of the ground state wavefunc- 
tion for TV particles, some theoretical work is needed to ex- 
tract experimentally relevant observables. As atomic den- 
sity profiles may be measured by absorption or even non- 
destructive imaging 12] , natural observables are functions 
of the positions of the particles. The simplest observable 
is the mean density of particles, p(x), obtained by an av- 
erage of the density profile over many (ideally infinitely 
many) experimental realizations. From the ground state 
iV-body wavefunction one however does not obtain any 
useful information on the mean density: Since the center 
of mass is fully delocalized, one gets a uniform distribution 
over the whole line. Experimental reality is very different, 
the soliton being obtained from an initially trapped Bose- 
Einstein condensate. When the trap is switched off to free 
the soliton, its center of mass is not in its ground state, it 
is in a localized and non-stationary state which depends 
on the experimental preparation procedure. A more real- 
istic assumption is thus to assume a iV-body wavefunction 
of the form 



${xi,...,x N ) =&(R)<f>(xi 



,x N ), 



where 



1 N 
N ^ 



(1) 



(2) 



is the center of mass position of the soliton and the center 
of mass wavefunction <P(R) is a priori unknown and de- 



Yvan Castin: Internal structure of a quantum soliton and classical excitations due to trap opening 



pends on the experimental details. The theoretical chal- 
lenge is thus to predict observables for a fixed position 
of the center of mass R. Experimentally relevant results 
are then obtained from these theoretical predictions by a 
further average over R with the probability distribution 
|^(i?)| 2 . One can even hope that the predictions for fixed 
R are measurable, if the center of mass position can be 
measured with a high enough accuracy for each individ- 
ual realization of the experiment. 

Turning back to the simple example of the mean den- 
sity, we see that the right concept is the mean density of 
particles p(x\R) for a fixed center of mass position R, as 
was already argued in [TTj within the general concept of 
(here translational) symmetry breaking. Remarkably an 
explicit expression of p(x\R) in terms of a sum of N — 1 
exponential terms may be obtained p~3|J14] . 

The next step beyond the mean density is the pair dis- 
tribution function of the particles p(x, y), or very similarly 
the static structure factor S(x,y), of the iV-particle soli- 
ton. Recently a general study of the dynamic structure 
factor was performed from the Bethe ansatz [§K], which 
includes the static structure factor S(x, y) and large N 
expansions as limiting cases. The goal of the present work 
is, in the spirit of the above physical discussion, to study 
the static structure factor S(x,y\R) for a fixed center of 
mass position R. This static structure factor gives access 
to correlations between the positions of the particles in- 
side the soliton, that is it gives information on the internal 
structure of the soliton, which goes beyond the usual mean 
field (or Gross-Pitaevskii) approximation, which neglects 
such correlations. 

As a guideline we imagine that, in an experiment, one 
wishes to access the variance of a one-body observable 
of the gas involving some function U(x) of the particle 
position x: 



w(R) 



■V 



E^ 



A' 



£^« 



(3) 



where the expectation value (. ■ -)r is taken over the in- 
ternal wavefunction <j){x\ , . . . , xn) of the quantum soliton 
for a fixed center of mass location R. 

The paper is organized as follows. In section [1] we re- 
call the basic facts about the ID Bose gas model and we 
give a general expression of the pair distribution function 
p(x,y\R) for fixed R that may be used to evaluate w(R) 
numerically for a moderate number of atoms and that will 
be the starting point of analytical calculations. In section 
[2] we present an analytical calculation of w(R) for an ar- 
bitrary atom number N, in the limit where the function 
U(x) is slowly varying over the spatial width of the soli- 
ton. In section [3] we present a large N expansion of the 
static structure factor S(x,y\R) for fixed R that we then 
use to calculate w(R) to leading order in N; in the same 
section, we integrate S(x, y\R) over R to see if we recover 
the results of [7] for S(x, y), and we also consider the case 
of a function U(x) much narrower than the soliton. In sec- 
tion 2] we evaluate the accuracy of the assumption (fT]) in 
an experiment: Assuming that the gas in the trap is in its 



ground state (at least for its internal variables), we cal- 
culate analytically (to leading order in N and in the trap 
curvature) the number of internal excitations of the soli- 
ton produced by the trap opening. We conclude in section 

El 



1 Model, basic definitions and general results 

1.1 Hamiltonian and ground state wavefunction 

We consider a set of TV > 2 spinless non-relativistic bosons 
of mass m moving on the one dimensional real line, in 
the absence of any trapping potential, with open bound- 
ary conditions, and binary interacting via an attractive 
contact potential of coupling constant g < 0. This corre- 
sponds to the Hamiltonian in first quantized form: 



H 



N 

E 



1m 



E 

l<i<j<A 



gS(xi 



(4) 



For this problem there is full separability of the center of 
mass motion and of the internal coordinates, the inter- 
nal wavefunction being independent of the center of mass 
state. As a consequence, the ground state wavefunction 
has it center of mass with zero momentum and depends 
only on the relative coordinates of the particles. Its exact 
expression is [5] 



[Xi,. 



,xn) = A/"exp 



m\g\ 
' 2h 2 



/ j \ x i x j\ 



l<i<j<N 



(5) 



where the normalization factor Af will be specified later. 
The corresponding eigenenergy is [9] 



Eo(N) 



mg 
24/J2 



(N - 1)N(N +1). 



(6) 



This ground state is the so-called quantum soliton with N 
particles. 



1.2 Mean density, pair distribution and static structure 
factor for a fixed center of mass position 

The crucial concept in the present work is the expectation 
value of a position dependent observable, that is of an ar- 
bitrary function 0(xi, . . . , xn) of the N particle positions, 
for a fixed value R of the center of mass position: 



(0)i 



/ dxi . . . dx N 5 I R — — y^ Xk J 
xO(xi, . . .,xn)\<P(xi, . . .,x N )\ 2 , 



(7) 



where the integration is over the whole space M. N . The 
soliton wavefunction shall then be normalized in such a 
way that the expectation value (1)r of the function O 
constant and equal to unity is also equal to unity. Using 
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the bosonic symmetry of the wavefunction, the integral 
to compute is then simply AH times the integral over the 
so-called fundamental domain of ordered positions 



D = (xi, . . . , xn) such that x\ < . . . < xn- 



(8) 



Over this domain, the wavefunction indeed has the sepa- 
rable expression 

</>( Xl ,...,x N ) =Afe-w£f= 1 ^-( JV + 1 )K (9) 
With the change of variables (fM| one then obtains 

-Jl x N-l 



WY 



m\g\ 



N 



(N-iy. 



i. 



(10) 



Setting O = ^ • 8{x — Xj) we obtain the mean density 
p(x\R) for a fixed center of mass position. This was first 
calculated in [T3l: 



p(x\R) 



N-2 

E 



AH ; 



N * t, (* 



(-l) fc (fc + l) 
2-k)\(N + k)\' 



-(k+l)\x-R\/Z 



where £ is the spatial width of the classical (that is mean- 
field) soliton, 

Z = -Tw- (12) 

m\g\N 

It is a function of \x — R\, a consequence of translational 
and parity invariance of the Hamiltonian and of <f>. One 
thus has p{x\R) — p(x — R\0) — p(R — x\Q). It is nor- 
malized in such a way that the integral over x over the 
whole space is equal to N. If the true physical state of the 
system is a product ([T]) of a localized center of mass wave- 
function <P(R) and of the quantum soliton wavefunction, 
the physical mean density is the average (or equivalently 
the convolution) 



p(x) = 



+ 00 



dR\<P{R)\ 2 p{x-R\0). 



(13) 



Setting O — J2i=£j $( x ~ x i)Hv ~ x j)> we obtain the 
pair distribution function p(x, y\R) for fixed center of mass 
position R. It is normalized as N(N — 1) for the double 
integration over x and y. From the translation invariance 
p(x,y\R) = p(x — R,y — R\Q) so that it is sufficient to 
calculate this pair distribution function for R = 0. In real 
space, we do not know an expression of p(x, y\0) as simple 
as (|llj) but we have found a simple expression for the 
Fourier transform 



p(q a ,q b \0)= / dxdyp(x,y\0)e-^ x +^. 

JR 2 

As detailed in the appendix [A~l one has 

r(N) 



(14) 



P(la,qb\0) = 



E 

l<j<k<N 



r(N + iQ) 

r(k-iQ)r{N + l + iQ-j) 

r(j)r(N + i-k) 



\r(j - 


-e a )r(N + l + iQ-e a -k) 


r{k- 


- e a )r(N + l + iQ-e a - j) 



£a <-> £fc 



(15) 



where r(z) is the Gamma function of complex argument z 
and we have introduced the dimensionless variables Q a .b = 
q a ,bi, Q = Qa + Qb and q = (Qb - Q a )/2. The quantities 
e a ,b solve the degree two equations: 



e a ,b(N + iQ- e a . b ) = iQ a ,bN. 



(16) 



In the large N limit, for fixed values of Q a and Qb-, it is 
convenient to take 



e a ,b = \{N + iQ)-[{iq±N/2f 



QaQb} 1/2 



(17) 



with a determination of the square root such that e a .b — 
iQa,b + QaQb/N + 0(1/N 2 ) for N ->■ +oo, e.g. with the 
line cut on the real negative axis. 

This gives the idea a posteriori to look for a similar 
expression for the Fourier transform p(q a \0) of p(x\0). As 
shown in the appendix[X]one obtains the simple expression 



(11) P(q a \0) 



r(N) 



r{N + iQ a ) 

y y r(j-iQ a )r(N + i + l Q a -j) ^ (18) 



r(j)r(N + i-j) 



Amusingly, this allows to express p{q a \0) in terms of the 
standard hypergeometric function evaluated in z = 1, 



p(q a \o) 



r(N)r(i-iQ a ) 
r(N - iQ a ) 

2 Fi(l-N,l 



iQ a ;l-N-iQ a ;l), (19) 



although we have not found this expression particularly 
useful. 

Let us come back to the original problem of calculating 
the variance w(R) of the one-body observable J^i U(xi), as 
defined in ([3]). One first expresses the squares as products 
of double sums over indices i and j. One can split the 
double sum over i and j in ([3]) in diagonal terms i = j 
and off-diagonal terms i ^ j, so that w(R) = Wding(R) + 
w fi (R) ■ The diagonal terms can be expressed in terms of 
the mean density, 



N 

E< 
t=i 



dxp(x\R)U{x) 2 - 



N 






dxp{x\R)U{x) 



(20) 



This can be directly evaluated from (p~Tj) . The off-diagonal 
terms depend on the pair distribution function, 

w oS {R) = Y J (U{x i )U(x j )) R 



dxdyU(x)U(y)p(x,y\R). (21) 
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Introducing the Fourier transform of U(x), 

U(q a )= / dxU(x)exp(-iqx), (22) 

Jr 

we obtain the Fourier space expression 
lOoff(-R) = 



This allows to express the Fourier transform of S(x, y) in 
terms of the Fourier transform of p(x, y|0), when one uses 



S{q a ,q b ) = 2ir5(q a + q b ) [N + p(q ai q b \0)} 



(29) 



dq a dq b - ~ 

-U{q a )U{q b ) 



(27T) 



xe ^+^ R p(-q a ,-q b \0) (23) 



From (|15[) we thus have an analytical expression of S(q a , q b ) 
in terms of a double sum. 

The relation (|28|) will also allow us, in the large N 
limit, to convert our large N expansion of S(x,y\R) into 
a large N expansion of S(x, y), see 33.21 



that can be directly evaluated using (fTSf . Simpler analytic 
expressions of w(R), either for slowly varying functions 
U(x) or in the large N limit, shall be given in the next 
sections. 

To conclude this subsection, we note that w(R) has a 
very simple expression in terms of the correlations con- 
tained in the fixed R static structure factor S(x, y\R), 

S(x,y\R) = (p(x)p(y)) R 

= S(x-y)p(x\R)+p(x,y\R) (24) 

where the operator giving the density is 

N 



2 Value of w(R) for a broad function U(x) 

One supposes in this section that U(x) varies slowly over 
the length scale £ of the quantum soliton, as defined in 
(|12| . Then one rewrites ((3|) using the translational invari- 
ance, 



w(R) 



N 



J2 U{xi + R) 



N 



) -{J2u(x, i + R))l (30) 



and one expands 



p(x) = ^<5(xi -x). 



(25) 



j=i 



We note that the double integral of S(x, y) over x and y 
is equal to N 2 . Using again the translational invariance, 
one obtains the illuminating expression 



w(R) 



dxdyU(x + R)U(y + R) 

x[S(x,y\0)-p(x\0)p(y\0)]. (26) 



U( Xi + R) = U(R) + XiU'{R) + -x 2 U"(R) + ... (31) 

The constant shift U(R) has no effect in w(R). The linear 
term has also an exactly vanishing contribution, since by 
definition ((J2i=i x i) n )o = for all integers n > 1. Setting 

N 



The essence of the mean-field approximation is to neglect 
correlations among the particles, so that S(x,y\0) would 
essentially be approximated by the uncorrelated product 
p(x\0)p(y\0). The above writing clearly reveals that w(R) 
is sensitive to correlations that are beyond the mean-field 
approximation. 



1.3 The usual static structure factor 

The usual static structure factor is the spatial correlation 
function of the operator p(x) giving the density, 



o 2 = ]>> 2 , 

we thus obtain the leading contribution 
w(R) 



i[[/"(i?)] 2 (Var0 2 ) 



o ' 



(32) 



(33) 



with (VarO 2 ) = (O 2 )o ~ (^2)0 is the variance of O2 for 
a center of mass position fixed at the origin of the coordi- 
nates. 

It turns out that an exact expression may be obtained 
for this variance, as detailed in the appendix [B] It is the 
sum of three contributions, 



(VarO 2 ) 



S(x,y) = (p{x)p(y)) 



(27) 



m\g\ 



[Si + S 2 + S 3 ] 



where the expectation value is taken literally in the ground 
state of the gas, thus assuming a perfectly delocalized cen- 
ter of mass wavefunction <P(R) = 1. Because of the trans- 
lational invariance, it is a function of x — y only. 

This usual structure factor is deduced from our fixed- R 
one by integration over the center of mass position 



with 



S(x,y)= f dRS{x,y\R)= f dRS{x - R,y - R\Q). 

Jr Jr 



N N N 

»=2 j=2 k=1 

N N 

i=2 3 =2 

N 



(28) 



S 3 = 6^2 B l 



(34) 

(35) 
(36) 
(37) 
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We have introduced the symmetric matrix 
1 1 

TD 

ij ~ N [N + 1 - xDhx(i,j)][max(i,j) - 1] ' 



(38) 



defined over the index range 2 < i,j < N. This holds 
whatever the value of the atom number N > 2. 

In the large N limit, the contribution Si is dominant, 
simply because it contains more terms, and its asymptotic 
expression is evaluated by replacing the sums by integrals, 



Si 



4 

7p 



d.r 



ln(l 



In a; 



1 

N~ 3 



- x 
8tt^ 
3 



16C(3) 



This leads to the estimate 



w(R) ~ A£ 4 [U"(R)Y 



2f_ 
3 



4C(3) 



(39) 



(40) 



for a slowly varying potential U{x) in the N ^> 1 limit. 



3 Large iV limit of w(R) for U(x) of any width 

In this section we give a large N expansion of the static 
structure factor for fixed center of mass position, which 
allows to get the asymptotic behavior of w{R) in the large 
N limit. When specialized to a quadratic potential U(x) 
the general result reproduces the large N broad potential 
result of the previous section. As a first test of the result, 
we integrate S(x, y\R) over the center of mass position R, 
to see if we recover the results of [7] for the usual static 
factor S(x,y). As a second test of the result, we get an 
approximate expression for w(R) in the case of a narrow 
function U(x), to lowest order in the width b of U(x), 
first from the large N expansion of S(x,y\R) and then 
from a more general reasoning not relying on a large TV 
expansion. 



3.1 Asymptotic expression of S(x,y\R) and of w(R) 

As we have seen in (f26]h w(R) is directly related to the 
deviation of the fixed R static structure factor of the soli- 
ton from the uncorrelated form p(x\R)p(y\R). It turns 
out that this deviation may be easily obtained from the 
Fourier space expressions (|15ll8p . simply by taking the 
large N limit of the r functions and by replacing the dis- 
crete sums over indices by integrals. This shows that our 
Fourier space representations are indeed useful. 

As detailed in the appendix [C] the large N expansion 
for a fixed value of £ (that is for a fixed value of Ng), gives 
the leading term 



5(a:,p|0)-p(a:|0)p(p|0) 



% - x) 



yji -xji-e 



x/i 



-Nd x d y 
-v/t 



[2cosh(x/2£)] 2 [2cosh(y/2£)p 



(41) 



where 9(x) is the Heaviside distribution that is equal to 
zero for x < and to one for x > 0. If one wishes to have 
an expansion of the static structure factor only, one has 
also to expand p(x\0) in powers of N. From the real space 
expression (fTTj) . expanding each factorial in the large N 
limit for a fixed summation index k, we obtain, setting 
X = x/£, a result in agreement with [13] : 



p(x\0) 



N 
J 



1 



1 d 2 

NdX 2 



1 



[2cosh(A/2)] 2 ' 



(42) 



where one may check the normalization condition N = 
J R dx p(x\0). This gives the expansion of S(x, y\0) up to 
order N, for a fixed value of £. 

From (|24[) and (|4"Tj) we can directly obtain, in the large 
N limit, the pair correlations between the positions of the 
particles for a fixed center of mass position: 



6p(x, y\R) = p(x, y\R) - p(x\R)p(y\R). 



(43) 



One can indeed show that the distributions generated in 
(T4"l"j) by the derivatives oi0(x — y) and 0(y — x) with respect 
to x and y exactly cancel with the Dirac term appearing 
in (|24| . As a consequence the expression for Sp(x,y\0) is 
deduced from the right hand side of (|4T|) simply by ex- 
changing the order of the 9 distributions and of the oper- 
ator d x d y . After an explicit calculation of the derivatives 
with respect to x and y we obtain 



fy(x,y\0) 



N 
16f 



[(2 + \X - Y\) sinh(A/2) sinh(y/2) 



+2 sinh(|A - Y\/2)] / [cosh(A/2) cosh(F/2)] 3 , (44) 

with X = x/£ and Y = y/£. A contour plot of this func- 
tion reveals that it has an interesting structure, in the form 
of two valleys separated by a crest on the x — y line, each 
valley being elongated in the direction parallel to x = y 
and containing two local minima separated by a saddle 
point. For clarity we only show a plot of 5p along the line 
y = —x, restricting to x > by parity, see the thick solid 
line in Fig[TJ the minimum of this line then corresponds 
in the full x — y plane to one of the aforementioned saddle 
points. In the same figure, we also give the value of Sp 
for finite values of A, obtained by calculating the Fourier 
transform of (|101|) and (I103|) numerically. This shows that 
the large N limit is well approached with moderately high 
values of N already. 

We come back to (|41|) to obtain the large N equivalent 
of w(R) at fixed £. By repeated integration by parts in the 
double integral of (|26)l . and using the fact that 



Y - X 



[2cosh(A/2)] 2 [2cosh(y/2)] ; 



= d x d Y - 



Y - X 



(e Y + l)(e~ x + I)' 

we find [15j for a function U(x) not rapidly increasing at 
| x | = co (that is not increasing faster than a power law): 



/-foo Z'+OO 

dX / dY 
-oo Jx 



U"{R + X£)U"{R + Yi) 



2 + Y - X 



(e Y + l)(e- x + 1) 



(46) 
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Fig. 1. Cut along the line y = —x of the function Sp(x,y\0) — 
p(x,y\Q) - p(x\0)p(y\0). Thick solid line: Large N limit (ji4)l . 
Dashed line: Numerical result for A^ = 10. Dot-dashed line: 
Numerical result for A^ = 20. Thin solid line: Numerical re- 
sult for N = 40. Sp is in units of N/£ 2 , and the coordinate 
x is in units of £. The inset shows numerical evidence for the 
convergence of 5p(0, 0|0)£ 2 towards a non-zero value ~ —0.031 
at large N, which is compatible with the analytical prediction 
6p(0,0\0)£ 2 /N — > 0. In the inset, symbols are numerical data, 
and the dotted line is a quadratic fit used to guide the eye. 

One immediately sees that a linearly varying poten- 
tial U(x) gives a vanishes contribution to w(R), which is 
also obvious from the definition ([3]), since ((X)i X *)")-R — 
(NR) n by construction. If U(x) is a quadratic function of 
x, U(x) = U"(R)x 2 /2, one should recover Eg. ([ID]) , which 
is indeed the case if 



7T 

T 



2C(3) 



+oo 



dX 



+ DG 



dY ■ 



Y - X 



x 



(e> 



l)(e-* 



We have checked that this identity holds |16j . 



(47) 



3.2 Application: large .V expansion of the usual static 
structure factor 

The usual static structure factor S(x, y) is obtained from 
S(x,y\R) by integration over R, see (j2"8")l . Assuming with- 
out loss of generality that y — 0, we see from (j4"Tj) and 
from (145)) that one has to integrate over R quantities of 
the form (d x d y F)(x — R,R), where the function F corre- 
sponds in a first stage to the expression in between square 
brackets in (|4Tj) and in a second stage to the function over 
which dxdy acts in the right hand side of (T45|) . From the 
differential relations, taking x and R as independent vari- 
ables, 



(d x F){x - R, -R) = — [F(x - R, -R)] 



(48) 



(d y F)(x-R,-R) 



dx dR 



[F(x-R,-R)], (49) 



we see that the integral over R will cancel the derivatives 
d/dR. Using 



+oo 



dR 



1 



X 



(e R + l)(e- R + x + I) e x -l : 



(50) 



one gets 



S{x,0) 



+ OC 



N d 2 
£ dX 2 



dRp{x\R)p(0\R) c 

d 2 X(X + 2) 



9(X) 



dX 2 



oX 



1 



X «-► -X 



(51) 



where X = x/£. One can show that the 9{X) distribution 
can be exchanged with the last d 2 /dX 2 if one wishes, but 
we shall not use this property here. From the large N ex- 
pansion (|42|) of p(x\R), and using the fact that the second 
order derivative of Xj [exp(X) — 1] is an even function, one 
finally obtains the large N expansion 



S(x,y = 

N d 2 
+ £ dX 2 



0) 



iV 2 d 2 



£ dX 2 



X 



,x 



-1 



6(X) 



X 2 



dX 2 e x - 1 



X 



-X 



(52) 



To compare this result with the ones of [7j given in Fourier 
space, one performs the Fourier transform of Eq.(62) and 
Eq.(63) of [7]. Then one finds that (|5^)) is consistent with 
[7J if one adds in Eq.(62) of [7J the subleading term of 
the large N expansion of the intermediate quantity S^(k) 
introduced in [7J. 

Furthermore, ([52]) can be shown to be equivalent to the 
more appealing writing, where the expected Dirac delta 
contribution is singled out, 



S(x,y = 0) = 



N 2 d 2 



£ dX 2 



X/2 



tanh(X/2) 



£ \ V 'dX 4 



X 2 /2 



tanh(X/2) 



X 



+ N6{x) 
-X 



(53) 



where we recall that X — x/(,. One can also easily check 
that this obeys the sum rule J R dx S(x, y — 0) = A^ 2 . 



3.3 Approximation of w(R) for a narrow barrier 

In this subsection, one assumes that U(x) and its deriva- 
tives are functions localized around the origin with a width 
much smaller than £. E.g. U(x) is a Gaussian centered in 
x = with a width -c £. One first considers the large N 
limit. We rewrite Eq. (|46|) as 



w{R) ~ 2N£ 



2 ' dx 



dy 



U"{x)U"{y)F 



x — R y — R 



with 



F(X 1 Y) = 



2 + Y - X 



( e Y + l)( e - x + l) 



(54) 



(55) 



Then we expand the factor containing F in powers of x/£ 
and y/£. If U(x) and its derivatives are rapidly decreasing 
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functions, one can show that the first non vanishing con- 
tribution to w(R) comes from the third order expansion, 
which results in 



w(R) ~p(R\0) 



+ OO 



dxU(x) 2 , 



(56) 



where we have used (|42l) to leading order in N to recognize 
the factor p(R\0). 

There is actually a faster way to obtain this result, and 
not restricted to the large N limit, from the exact writing 



dx p{x\R)U (x) 2 + I dxdy 

Jr 2 

[p(x, y\R) - p(x\R)p(y\R)]U(x)U(y). (57) 



w(R) 



Then one sees that the first term in (|57jl is first order in 
the width b of U(x), since the integration range over x has 
a width 0, whereas the second term is second order in b, 
since it involves a double integral over a range of diameter 
~ b. Then to first order in b one recovers (1561) . In the large 
N limit, for a fixed £, both terms of (|57|) scale in the same 
way with N, that is linearly with TV, so that the order of 
the limits N — » +00 and b — > (for fixed £) is not crucial. 



4 Number of internal excitations created by 
the trap opening 

In this work, we have assumed up to now that a pure 
soliton is produced in the experiment, corresponding to 
the state ([T]): The center of mass may be in an arbitrary 
excited state but the internal variables of the gas are in 
their ground state. In this section, we revisit this assump- 
tion taking into account experimental constraints. 

In a real experiment the ultracold gas is prepared in 
a trap; it is not free along x, and each atom is subject to 
the harmonic confining potential 



W(x) = —muj x . 



(58) 



We recall that the center of mass motion and the internal 
variables remain separable in a harmonic trap. In order 
for the trapped gas to be close to the free space limit, 
the oscillation frequency u is adjusted to have Hu> <C \po\, 
where 

h 2 . . 

am€ 4 

is the mean-field approximation for the chemical potential 
of the free-space soliton. We assume that the trapped gas 
is cooled down to a temperature T low enough to have 
ksT <C \po\- The gas is then a pure soliton. On the con- 
trary we do not assume the much more stringent condition 
ksT -c hu, so that the center of mass of the gas may still 
be in an excited state. 

Eventually the trapping potential along x will be swit- 
ched off, to obtain the ideal conditions of the Hamiltonian 
(|4| . This trap opening will create some internal excitations 



of the gas, that is the untrapped gas will not be a pure 
soliton but will contain a mean number N cxc of internal 
excitations. 

Since no Bethe ansatz solution exists in a trap, the 
modest goal here is to calculate N exc to leading order in N 
in the large N limit, for a fixed value of \po\/(huj) 3> 1 |17j . 
In this limit, a ID classical field treatment is sufficient, 
with a Hamiltonian 



H 



dx 



-^* d & + |V>*V + \muJ 2 (t)x 2 r^ 



(60) 
We have slightly generalized (|58|) to treat the case of a 
time dependent trap: 



w 2 {t) = u 2 X {t) 



(61) 



where the function x(i), going from unity for t = to 
zero for t — > +00 describes the switch-off procedure of the 
trap. The atom number is fixed so that the norm squared 
of the classical field is also fixed: 



dx\tp{x)\ 2 =N. 



(62) 



This classical field problem can be solved on a computer, 
which will allow a test of our predictions for N cxc /N. 

The problem can be further simplified in the limit 
Hlo/\po\ —> 0, where one expects that the field ip{x) will 
remain "close" (up to a phase factor) to the one describing 
a pure soliton [IS], il>o(x) — A rl / 2 ^ (x), with 



<t>o{x) 



1 



1 



2^/2 cosh[x/(2£)] ' 



(63) 



We then use the number conserving Bogoliubov formalism 
of [TW2T)] , downgraded to a classical field problem (simply 
replacing commutators with Poisson brackets) , as was al- 
ready done in [21]. One splits the field by projection along 
the mode 0o and orthogonally to it: 



rp(x) = a (j> (x) +tp±(x) 



(64) 



where ao is the component of the field on the mode </>o 
and the field ip±(x) is orthogonal to that mode. The idea 
is to treat ip±_ (x) as a small perturbation. The strength of 
the number conserving approach is to eliminate the am- 
plitude ao in a systematic way, using the modulus-phase 
representation 

ao = \ao\e 10 . (65) 

The phase is eliminated by a redefinition of the trans- 
verse field: 

A(x) = e-*Vj.(aO. (66) 

The modulus |ao| is expressed in terms of A using the 
condition of a fixed atom number (|6"2"|) . 

In the absence of trapping potential, one keeps terms 
up to quadratic in A in the Hamiltonian, and the resulting 
quadratic form can be written in normal form as [22j 



Ti-o — £0 



P 2 



2mN 



dk 

2^ 



ekbtbk 



(67) 
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where £q = N/j,q/3 is the ground state of the classical field where (differently from section[3]) we have set X = x/(2£) 



model. The variable P of the field represents the total 
momentum of the field, written to first order in A: 



P 



-A 1 / 2 



dx<j>' Q (x)[A*{x)-A(x)}. 



(68) 



and K = 2k£. The modes for k < are deduced from the 
relations Uk(x) — U- k (—x) and Vk(x) — V- k {—x). 

In presence of the trap, there is an extra contribution 
to the Hamiltonian containing the trapping energy, 



The occurrence of the term P 2 / (2mN) represents phys- 
ically the fact that the center of mass motion is decou- 
pled; more formally, it corresponds to the fact that <j>q (x) 
"spontaneously" breaks the translational symmetry of the 
Hamiltonian, which leads to the occurrence of a Goldstone 



W= -mu 2 (t) 
2 v ' 



dxx ip*(x)ip(x). 



(78) 



We approximate it to leading non-trivial order in A, that is 
to first order, injecting the splitting Ij64|) . The term in a^ao 
deviates from N by 0(A 2 ) so it gives to first order only 
a constant contribution, with an integral over x that can 



mode ggmm ■ Th e field variable canonically conjugated be calculated exac% if ' neC essary. The crossed terms a* A 



to P corresponds to the center of mass position of the field 
written up to first order in A 



Q = N- 1 / 2 f dxx<j> Q {x)[A{x) +A*(x)]. 



(69) 



and complex conjugate are kept, with ao approximated by 
iV 1 / 2 e l9 , whereas the quadratic terms A* A are neglected. 
Finally, we replace the field A{x) by its modal expansion 



Apart from this Goldstone mode, the other eigenmodes 
behave as a continuum of decoupled harmonic oscillators, 
with normal (complex) variables b k and eigenenergies 



f-k = \m\ 



h 2 k 2 

2m 



(70) 



A(x) = -N 1 / 2 ^{x)Q + 



We thus keep 



HN 1 / 2 



X(f> (x)P 



dk 

-[u k (x)b k +v* k (x)bl]. (79) 



W 



which correspond to internal (and thus gapped) excita- 
tions of the gas: An elementary excitation physically takes 
the form of a free particle coming from infinity with a 
wavevector k and scattering on a soliton with N — 1 par- 
ticles. The field variable b k has the expression 



-Nmu) 2 (t) 
2 v ' 



dxx 2 cj) 2 t (x)+x(t) 



with coefficients 
Ik = N 1/2 -m^ 



^hkbk+l* k b* k }, 
(80) 



dxx 2 (j) (x)[u k (x) +v k (x)}. (81) 



b k 



dx [ut(x)A(x)-v* k (x)A*(x)} 



(71) 



With the residues method, the resulting integral can be 
evaluated analytically. We give for convenience the ratio 
to the mode eigenenergy: 



where the Bogoliubov modes of the number conserving 
theory are expressed as follows in Dirac's notation (see 
[20 §V.A) here for a real function 0q: 



\u k ) 

\Vk) 



Q\U k ) 
Q\V k ) 



(72) 
(73) 



Tfe _ 
where K 



IM) 



w(N£) 1 / 2 /2 



1 



(K - i) 2 (l + K 2 ) cosh{K-K/2) 



(82) 



where Q = 1— |<fo)(0o| projects orthogonally to \<fio), and 
the Uk, V k are the eigenmodes of the usual Bogoliubov-de 
Gennes equations 



2fc£ as above. 

When the perturbation W is added to the unperturbed 
Hamiltonian (|67p , the equilibrium value of the b k 's mini- 
mizing the resulting energy is 



MO) = -— ■ 

Cfe 



(83) 



ekU k (x) = 



-e k V k (x) = 



h 2 
2m" 



+gN<f> 2 (x)V k (x 

h 2 



2gN\<f> (x)\ 2 ~ »o 



U k (x) 



-—d 2 x + 2gN\Mx)\-Vo 
-gN4>* 2 {x)U k (x). 



V k (x) 



It turns out that these modes are known exactly 
k > one has 



(74) 
i 

(75) 
. For 



At later times, the trap is opened according to the switch- 
off function x(i). The Hamiltonian equations of motion 
then give 

ih di bk ^ = tkhk ^ + xW^' 
to be solved with the initial conditions 

el 



(84) 



b k {t) = b k (0) 



X {t)-e- 1 ^ 1 I dre 1 ' 
o 



■d X (r) 



dr 



(85) 



U k (x) = e iKX 

V k (x) = 



1 + (K 2 - 1) cosh 2 X 



(K 



iKX 



{K -i) 2 cosh z X 



where Lu k — e k /h is the mode frequency. This allows to 
calculate the number of excitations after the trap was 
2% K sinh X cosh X switched off, using x(£) — > for t — > +oc: 

i) 2 cosh X f dk f dk 

N m = Um / ^bt(t)b k (t)= / -|M0)| 2 /H) 

{ll) (86) 
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where I{D) is a spectral density of the switch-off proce- 
dure at frequency f2: 



I{Q) 



+ OC 



dr e 



-in 



T dx(r) 



dr 



(87) 



We recall that (|86|) is valid up to leading order in N (be- 
cause of the classical field model) and to leading order in 
uj 2 (because of the perturbative treatment of the devia- 
tions of the field from the free space soliton). 

The case producing the maximal number of excitations 
for a monotonic x(i) corresponds to a sudden trap switch- 
off, where I{£2) = 1 at all frequencies. Using (|82|) . we 
find that the resulting integral can be evaluated with the 
residues method, so that 



N, 



sud = CN 






with 



16 



dK- 



(1 

7T 2 (7T 2 + 25) 

3840 



K 2 )^ cosh 2 (Kit/ 2) 
C(5) + 5C(3) 



128 



= 0.1446785. 



(89) 



This result is encouraging since a moderately small value 
fku — | /xo |/10 already leads to a number of excitations 
relative to the total atom number at the 10~ 5 level. As 
shown in Figj2] the analytical prediction (f88|) is in good 
agreement with the number of excitations deduced from a 
numerical solution of the Gross-Pitaevskii equation, pro- 
vided that |/j,o| 3> fuo. 

The number of excitations can be reduced by switching 
off the trap more slowly. E.g. a linear ramping 



X (t) = (l-\t)6(l-\t) 



leads to 



I{Q) 



l [n/{2\)\. 



(90) 



(91) 



Due to the presence of a gap |/xo | in the internal excitation 
spectrum of the gas, one gets the upper bound on the 
number of excitations 



f2h\ 
V Mo 



iVC 



(92) 



that is one gains quadratically with the switch-off time 
when it becomes longer than the internal soliton time 

h/\n Q \. 




Fig. 2. For an initially harmonically trapped classical soliton, 
number of excitations N^c produced by a sudden opening of 
the trap, divided by the number of particles N, and given as a 
function offtw/|/xo |, in log- log scale. Solid line: Analytical result 
(|88[) obtained in the limit |/xo| ~S> fou). Symbols: Result deduced 
from a numerical solution of the Gross-Pitaevskii equation in 
the trap. Here ui is the oscillation frequency of the particles in 
the trap, and /io = —mg 2 N 2 /8h 2 is the chemical potential that 
the classical soliton would have in the absence of trapping po- 
tential, m being the particle mass and g the coupling constant 
describing the interactions in ID. 



particles in a quantum soliton for a fixed position of the 
center of mass of the soliton. In particular, we have ob- 
tained analytically the large N limit expression of these 
pair correlations, see (|4Tj) and (|44|) . By integrating (J4TJ) 
over the center of mass position, we obtain a large TV ex- 
pansion of the static structure factor for a fully delocalizcd 
center of mass position, a quantity already studied in [6l 
[7] . On an experimental point of view, our predictions can 
be tested by measuring the positions of the particles in a 
very broad quantum soliton, prepared with weak attrac- 
tive interactions and a relatively small atom number |27j . 

The second problem was studied in the classical field 
model. The number of internal excitations of the gas cre- 
ated by the trap opening from an initial pure soliton was 
calculated in the limit where the soliton size is smaller 
than the size of the harmonic single particle ground state, 
see (|88l) for a sudden trap opening. This also can be seen 
experimentally by detecting atoms flying away from the 
remaining soliton core after the trap opening. A possible 
extension of this calculation is to include quantum fluctu- 
ations of the field. 

We acknowledge useful discussions with L. Khayko- 
vich, C. Weiss, A. Sinatra, M. Olshanii. Our group is a 
member of IFRAF. 



5 Conclusion 

Inspired by recent observations of matter wave bright soli- 
tons in atomic gases, we have considered here two prob- 
lems that may be relevant for experiments. 

The first problem is strictly beyond mean field: It cor 



A Calculation of the Fourier transforms of 

p(x\0) and p(x,y\0) 

We start with the definition of the mean density for a fixed 



responds to the pair correlations between positions of the center of mass position R, taking here R = without loss 
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of generality: From (|7|) one has 

p(x\0) 



/ dxi . . . dx N 6 y^ x k /N 

Jrn Vfci / 



A 



E 6 ( x j - x ) 

j=x 

x\(j ) (x 1 ,...,x N )\ 2 . (93) 



where we used the fact that the product of all (3 k (for k 
from 2 to JV) is equal to [(JV — l)!] 2 . Replacing the a k .j 
by their expression and using the identity 



N-x 



H(k+z)= 



k= 3 



r(j + z) 



Vj e i, . 



,N 



(102) 



Using the bosonic exchange symmetry we can restrict the 
integral to the fundamental domain D of ©, including a 
factor N\. For simplicity we take in this appendix h 2 /(m\g\) 
as the unit of length. With the change of variables, of Ja- 
cobian equal to unity, 



J2u k , l<j<N, 



deduced from the basic property r(z + 1) = z r(z) of the 
Gamma function, with the convention that an "empty" 
product is equal to unity, one gets fT8"]). 

We now turn to the pair distribution function for a cen- 
ter of mass position fixed in R = 0. Following the same 
steps as for the mean density, we obtain the Fourier trans- 
form 



(94) 



N 



fe=i 



K<la,qb\0) 



the condition to be in D is simply that all the U2 , ■ ■ ■ , ujv 
are positive, and u\ can vary in the whole real space R. 
Setting 



e n 

X<j=£k<N n=2 



Pn 



fln + iq a a n ,j + iqbOi n ,k 



(103) 



A' 



(3 k = YfiJ -(N + l)] = (N + l-k)(k- 1), (95) 

j=k 



we obtain from ©: 

P (x\o) = i(N-iy.fJ2[ du J' 



We can restrict to a summation over j < k, the contri- 
bution for j > k being deduced by exchanging q a and 
qt,. In the product over n, three ranges have then to be 
considered, (i) the first range n < j, (ii) the mid-range 
j + 1 < n < k, and (iii) the last range k + 1 < n < N. 
The first range and the last range contributions can be 
expressed in terms of the Gamma function as in (|102[) . 
The mid-range contribution is equal to the product 



duo . . .du 



N 



(N - n')n' 



3=X 



k-X 

mid = -W (N - n')ri - iQJN - n') + iQ b ri 

n'=j ' 



(104) 



s(f: N+ ^ k u k )s(x-J2^)e-^^^. (96) 



vfc=l 



k=X 



One can calculate the integral over u\: The first delta 
factor in |96|) . the one ensuring that R = 0, imposes a 

value 

"i = - S s? Uk - ( 9? ) 

fe=2 



where we have reindexed the product setting n' = n — 1 
and we have defined Q a — q a /N and Qb = qb/N. The last 
step is to consider the denominator in each factor of (|104|) 
as a polynomial of degree 2 in n': Its roots are e a defined 
in (HI EJ) and iV + i(Q a + Q 6 ) - e a . This leads to 



fc-i 



N 



^inid — 



11 M 



(JV - «')«' 



n =3 



(»' - e a )[JV + i'(Q a + Qb) - e„ - n'] 



7T ( 105 ) 



When one reports this value of Mi in the argument of the 
second delta factor in (|96p , one obtains a remaining factor 



• A' 



•K* ~ J2k=2 a K 3 u k), with 
ifc-1 



which can now be expressed as a ratio of products of 
Gamma functions. Then one gets ([15)) . 



"fe,.7 = 



JV 

7V + 1 -k 



N 



for fc < j 

for fc > j. 



(98) 
(99) 



Taking the Fourier transform of this remaining delta fac- 
tor, according to p(q a \0) — J R dx e~ tqaX p(x\0) , gives 



N 
p(q a |0) = [(JV-l)!] 2 ^/ du 



. du 



N 



3=1' 



JV JV 

En 

j = l fe=2 



/4 



/3fe +iq a a k , 



(100) 
(101) 



B Calculation of the variance of 2 

We explain how to calculate exactly the moments (6*2)0 
of the quantity O2 defined in ([32l in the JV-body internal 
ground state for a fixed center of mass position R = 0. We 
take h 2 /m\g\ as unit of length and we use the transforma- 
tions exposed at the beginning of appendixfA] Considering 
the change of variable (|9~4"|) , we rewrite O2 as 

JV 

C 2 =^[(^-Ml)+Ml] 2 
i=X 

I JV \ JV 

= -JVw 2 + 2w a I ^Xi I -(-^(^-u!) 2 . (106) 



i=l 
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The first sum in the right hand side of (|106p will have 
a vanishing contribution, since the expectation value is 
taken for a zero center of mass position. Replacing Xi — U\ 
by its expression in terms of U2, ■ ■ ■ , mjv, and using the fact 
that the value of u\ is fixed to ((57)) we see that one may 
effectively replace O2 by the quantity 

»=2 \fe=2 / \fc=2 / 

N 

— 5Z AijUiUj (107) 

with the symmetric matrix 

Aij = i[JV + l-max(t,j)][min(i l i)-l]. (108) 

We have thus reduced the problem to the calculation of 
the integrals 



(0 2 l )o = {(N-l)lY I du 2 ...du N 



k. It thus makes sense to split p(q , a |0) / 5(<7f ) |0) into a off- 
diagonal part (j ^ k), that we collect with the double 
sum in p(q a , q b |0), and a diagonal part 



Diag 



N 



n 



j=l ?7=a,b 



r(N) 



r{N + iQ v ) 

r(N+l + iQ v - j)r(j - iQr,) 

r(N + i-j)r(j) 



(112) 



where Q a , b are defined above P^|) . 

The idea to obtain the large N limit is simply to re- 
place the discrete sums by integrals. To this end, one has 
to calculate the large N limit of each term of the sums, for 
fixed values of y a = j /N and y b = k/N. The expansion of 
the Gamma functions is conveniently performed using 



r{z + b) 



{a-b)(a + b-l) 

Yz 



+ 0(l/z 2 ) 

(113) 

where the real quantity z tends to +00, and the fixed 
quantities a and b may be complex. One also uses the 
large N expansion of the quantity e a : 



(114) 



JY 



y^ A^muj 1 1 — '■ ~ 



-Y.k = 2PkUk 



(109) 



\hj=2 



£a = iQa H Tj h • • ■ 

For the diagonal part, only the leading term of (|113p is 
useful. Replacing ^ by N J dy a leads to 



with (3k defined in (|95|) . These integrals may be calculated 
by interpreting them as Gaussian averages, introducing 
the auxiliary complex random variables otk, 2 < fc < N: 
These variables ctk are statistically independent and each 
one has a Gaussian probability distribution ex e~l aA: l @ h . 
Each ut then corresponds to |ckfc| 2 , so that 



Diag : 



N I dy a e l 
/o 



Mn- 



= N dX a 



aiQX a 



2 cosh ^) 2 ' 



A 



(O 2 l )o = {([ J2 A > 



1 i2 1 i2 



))• 



(110) 



(115) 
where Q = Q a + Qb and the integral was transformed with 
the change of variable X a — In ,^° , to acquire the form 
of a Fourier transform. Note that the resulting integral 
can be calculated exactly, giving 7rQ/sinh(7r(5), but this 
is not useful here. 

For the double sum over j < k, one has to include 
the 1/z term in the expansion (|113|) . to obtain a non-zero 
result: 



Here ((. . .)) denotes the Gaussian average of the ct^'s and 
can be calculated using Wick's theorem. The calculations 
are a bit lengthy for n — 2 since Wick's theorem has to 
be applied to a product of 8 variables. Clearly the matrix 
Aij/ (Pi/3 j) appears, which is the matrix Bij of (|38|) . We 



Off-Diag ~ NQ a Q b / dy a 
Jo 
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dy b e l 



. In- 



V Qfcln ^f 



finally get 



for(Var0 2 ) = (Ol)o-((0 2 )o)^ 
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Uu 
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+ Qa 
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(116) 



C Large N asymptotics of 

p(x,y\0)-p(x\0)p(y\0) 

To perform the large N expansion we shall use the Fourier 
space expressions (I15ll8p . We thus define 



The change of variables X„ = In . Va and Xh = In . Vh 
gives to the off-diagonal contribution the form of a Fourier 
transform: 



g— iQa.X a e — iQbXt, 



Sp{Qa,qb\0) = p(q a ,qb\0) - p(q a \0)p{q b \0). 



(Ill) 



Off-Diag ~ NQ a Q b / dX a dX b 

JR 2 

x 6(X b ~ X a ) [X b -X a - e x « - 



(2cosh^) 2 (2cosh^) 2 



-Xi 



y b - 

(117) 



The pair distribution function involves a double sum over 
j ^ k. Since p(q\0) is a simple sum, the last term in (II 1 ip 
is a double sum over j and fc, without the restriction j =/= 



The leading term of Sp(q a ,q b \0) for N — > +00 for a 
fixed £ is the sum of (j!15[) and (|117[) . The Fourier trans- 
form with respect to q a ^ b is straightforward: The factors 



12 



Yvan Castin: Internal structure of a quantum soliton and classical excitations due to trap opening 



Qa.b act as derivatives, and the remaining bits have al- 
ready a Fourier form. The Fourier transform of the diag- 
onal contribution gives a contribution involving a factor 
8{x a —Xf,), which exactly cancels with the first term in the 
right-hand side of ([M)l . at the considered order in N, see 
(El. We obtain BIT). 
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